import sys,os 

whatshap = "/scratch/project_mnt/S0030/jingping/mamba/envs/whatshap/bin/whatshap"
minimap2 = "/scratch/project/qaafi-cnafs/jingping/software/miniconda3/bin/minimap2"
ref = "/scratch/project_mnt/S0030/jingping/02.Pineapple_genome_project/14.switch/callsnp/ref/Hap1.fa"
bwa = "/scratch/project/qaafi-cnafs/jingping/software/miniconda3/bin/bwa"
sambamba = "/scratch/project/qaafi-cnafs/jingping/software/miniconda3/bin/sambamba"
fq1 = sys.argv[1]
fq2 = sys.argv[2]

cmd = f"{bwa} mem -t 80 -M -R '@RG\\tID:stLFR001\\tSM:stLFR001\\tPL:ILLUMINA' {ref} {fq1} {fq2} | samtools view -bhS - |samtools sort -@ 80 -o stLFR001.bam\n" 
cmd += f"{sambamba} markdup -r stLFR001.bamwhic stLFR001.rmdup.bam\n"
cmd += f"gatk --java-options \"-Xmx200g\" HaplotypeCaller --tmp-dir tmp/ -R {ref} -I stLFR001.rmdup.bam -O stLFR001.g.vcf\n"
print(cmd)
#cmd += f"{minimap2} -t 80 -ax map-pb --secondary=no {ref} pacbio.subreads.fasta.gz |samtools view -bt {ref}.fai - |samtools sort -@ 80 -o pb.sorted.bam -\n"
#cmd += f"{whatshap} phase --ignore-read-groups -o Chr01HA.phased.vcf --reference={ref} --chromosome Chr01HA TGY.illuminaWGS.vcf pb.sorted.bam"